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I.  INTRODUCTION 


The  elliptic  coverage  function,  ELP,  defines  a  probability  function  P(r,  h,  k,  u,  v).  ELP 
gives  the  probability  of  a  shot  falling,  under  a  bivariate  uncorrelated  normal  distribution 
with  mean  zero  and  standard  deviations  u,  v,  in  a  circle,  T,  with  radius  r  and  centered  at 
(h,  k)  of  the  xy  plane.  This  probability  is  given  by 


P(r,  h,  k,  u,  v)  = 


fh+r  />k+  yj  r1 2  —  (x—  h)2 


2tT\\  V  Jh_r  Jl;--y/r2-(x-h)2 


exp  <  — 


,xN 


©2  +  (b2l  Mydx.  (i) 


2  L  u' 


v 


Introducing  dimensionless  variables,  let  x  =  (\/2u)  s  and  y  =  (\/2  v)  t,  then  (1)  becomes 


1  /-H+R  /•K+(u/v)v/R2-(s-H)2 

P(R,  H,  K,  u/v)  =  —  /  exp(— s2)  /  _  exp(— t2)dtds  (2) 

71  J  H-R  J  K-(u/v)^/R2-(s-H)2 

with  R  =  r/(y/2u),  H  =  h/(\/2u),  K  =  k/(\/2v).  Finally,  (2)  can  be  written  as 

i  /“H+R 

P(R,H,K,u/v)  =  — =  /  Fl(s)  ds,  (3) 

2  v71-  </h-r 

where 

Fl(s)  =  exp(— s2)  aerf  [K,  (u/v)a/R2  —  (s  —  H)2  ],  (4) 

,  o  /■a+b 

an  aerf(a,  b)  =  — =  /  exp(— z2)dz.  (5) 

V71  J  a— b 

Weapon  target  studies  ([3],  [9],  [10])  often  need  the  inverse  problem  solution  of  finding 
r  given  P,  h,  k,  u,  v.  An  objective  of  this  report  is  to  describe  an  algorithm,  INELP,  to 
find  r,  and  to  discuss  the  associated  Fortran  77  computer  program  INVELP  and  some  of  its 
supporting  routines. 

The  Fortran  source  file  IELP.FOR  contains  38  double-precision  routines,  which  are  item¬ 
ized  in  the  next  section.  They  include  the  main  subroutine  INVELP  and  all  its  required 
supporting  routines,  some  developed  by  the  author  and  the  rest  taken  from  [8].1  For  exam¬ 
ple,  the  computation  of  P  from  (3)  is  required  throughout.  It  is  carried  out  by  the  integration 
subroutine  DQXGS  of  [8,  p.  531].  2 

Four  appendices  are  included.  The  procedure  used  to  reduce,  whenever  possible,  the 
integration  interval  in  (3)  is  given  in  Appendix  A.  Appendices  B  and  C  contain  discussions 
for  obtaining  first  estimates  to  r.  Appendix  D  contains  a  table  for  r,  which  is  given  to  six 
significant  digits,  similar  to  the  five-digit  inverse  table  of  of  [2], 


1The  software  development  was  carried  out  using  an  IBM  PC  with  Lahey  Fortran  [6]. 

2  Ordinarily  PKILL  of  [8,  p.  127]  would  be  used  because  of  its  efficiency,  but  it  was  found  that  the  accuracy  set  for  determining 
r  could  not  be  obtained  with  PKILL  in  some  extreme  cases. 
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II.  ANALYSIS  TO  DETERMINE  r 

In  this  section,  the  analysis  used  in  the  algorithm  INELP  to  determine  the  radius  r  =  rQ 
of  the  target  circle  T,  given  P  =  PQ,  h,  k,  u,  v,  is  described  in  general  terms.  Greater  detail 
can  be  obtained  from  the  appendices  and  the  code.  To  simplify  notation,  P(r)  or  P  will  be 
used  to  denote  P(R,  H,  K,  u/v)  or  P(r,  h,  k,  u,  v),  with 

P(r)  =  P(r)  -  P.  =  P  -  P„. 

An  algorithm  for  the  numerical  evaluation  of  P  from  (3)  is  needed.  This  computation,  as 
noted  in  the  Introduction,  is  achieved  by  the  routine  DQXGS  of  [8],  which  uses  an  adap¬ 
tive  integration  scheme.  It  requires  routines  for  the  numerical  evaluation  of  Fl(s)  and  the 
aerf  function  as  noted  by  (4)  and  (5).  The  double-precision  functions  FI3  and  DAERF 
[8,  p.  51]  serve  this  purpose.  In  addition,  it  has  been  found  by  extensive  experimentation 
that  (3)  can  be  evaluated  numerically  more  quickly  if  u  <  v.  Consequently,  if  v  <  u,  then  u 
and  v  are  interchanged  as  well  as  h  and  k,  which  from  (1)  or  (3)  amounts  to  interchanging 
the  order  of  integration  leaving  P  unchanged. 

The  routine  DQ,  which  calls  DQXGS,  first  attempts  to  reduce  the  interval  of  integration 
[H  —  R,  H  +  R],  because  often  exp(— s2)  becomes  very  small  at  one  or  both  ends  of  the 
integration  interval  and/or  since  the  aerf  function  of  the  integrand  is  zero  at  H  —  R  and 
H  +  R.  See  Appendix  A. 

The  underlying  idea  for  finding  rQ  is  to  trap  it  between  a  smaller  value  Amin  and  a  larger 


value  Arnax,  so  that 

Amin  <  rQ  <  Amax.  (6) 

The  first  estimate  iq  for  rQ  is  obtained  from  (see  [2,  pp.  10-15]) 

D1  =  V  h2  +  k2  (7) 

D  =  D1  —  7  u  (8) 

ri  =  max[D,  (h v  +  k u  —  7  \/2  uv)/V u2  +  v2  ,  0].  (9) 


The  second  estimate  for  rQ,  r2  =  Rg,  is  obtained  from  Grubbs’  approximation  for  r0,  using 
the  routine  GRUB  [5].  This  estimate  is  often  very  good,  but  it  can  also  at  times  be  poor. 
Appendix  B  contains  Grubbs’  equations  taken  from  [5]. 

The  third  estimate  for  rQ,  r3  =  Rsq,  is  obtained  from  an  expression  for  the  probability, 
Psq,  of  a  shot  falling  under  a  bivariate  uncorrelated  normal  distribution  with  mean  zero 
and  standard  deviations  u  and  v  in  a  square,  S,  that  circumscribes  T  with  sides  parallel  to 
the  coordinate  axes  and  of  length  2  Rsq.  The  Newton- Raphson  (N-R)  procedure  [7,  p.  119] 
is  used  to  find  Rsq  with  Psq  =  PQ.  Consequently,  Rsq  <  A  although  truncation  and/or 

3  If  no  reference  is  given,  that  routine  was  developed  specifically  for  IN  YELP. 
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rounding  error  may  cause  this  inequality  to  be  violated.  The  (N-R)  analysis  is  given  in 
Appendix  C.  The  involved  subroutine  is  PSQR. 

At  this  point,  for  i,  j  =  1,  2,  3, 

Amin  =  max(si,  s2,  s3 ),  where  S;  =  q  if  P(ri)<0;  s;  =  0  otherwise,  (10) 

Arnax  =  min(ti,  t2,  t3  ),  where  tj  =  rj  if  P (rj)  >  0 ;  tj  =  1050  otherwise.  (11) 


For  example,  if  P(ri)  <  P(Rg)  <  P(Rsq)  <  PQ,  then  Rsq  is  chosen  as  a  starting  value  for 
Amin.  On  the  other  hand,  if  P(Rg)  >  P0,  then  Amin  =  Rsq,  and  Arnax  =  Rg  are  taken  as 
starting  values. 

It  is  also  easy  to  see  that  S  is  circumscribed  by  T  if  the  radius  of  T  is  taken  as  r  =  \/2Rsq. 
In  this  case,  if  \/2Rsq  <  Arnax,  then  Arnax  is  set  to  \/2Rsq. 

Another  estimate  is  made  to  obtain/improve  Arnax  by  the  result,  derived  in  [2,  p.  15], 


that 


ro  <  I'm 


D1  +  max(u,  v)  a/— 21og(l  —  PG) ,  if  PQ  >  10  10, 

D1  +  max(u,  v)  \J  2  PQ  ,  if  0  <  PQ  <  10~10. 


(12) 


Therefore,  if  rm  <  Arnax,  then  Arnax  =  rm.  At  this  stage,  values  for  Amin  and  Arnax 
have  been  established,  where  it  is  understood  that  associated  values  of  Pmin  =  P(Amin) 
and  Pmax  =  P(Amax)  have  also  been  determined  using  DQ.  In  fact,  at  any  stage  when 
a  new  estimate  r  for  rQ  is  found,  P(r)  is  computed.  In  addition,  if  at  any  stage,  with 
EPS3  =  10-10,  W1  =  1/10  if  P0  >  .99999,  else  W1  =  1, 

|P(r)  -  P0|  <  WE  =  EPS3  W1  PG,  (13) 

or,  with  WEI  =  10“14, 

Arnax  —  Amin  <  WEI  [Arnax  +  Amin]/2,  (14) 

then  an  acceptable  value  for  rQhas  been  found.  INVELP  exits,  with  outputs  (see  Call  line 
of  INVELP,  page  5)  R  for  rQ  and  P  for  P  =  P(r)  —  PQ. 

If,  at  this  stage,  neither  (13)  nor  (14)  holds,  which  is  generally  the  case,  then  the  interval 
Arnax  —  Amin  is  reduced  further  by  the  following:  Let 


Da  =  Arnax  —  Amin.  (15) 

Then 
NJ  =  7 

DaNJ  =  Da/NJ 
IF  (Pmax  >  —Pmin)  THEN 
DO  25  J  =  1,  NJ 

A  =  Amin  +  J  *  DaNJ 

CALL  DQ(A,  P,  I)  ! Routine  to  compute  P(A) 

IF  (P  >  0)  THEN  !I  =  No.  of  calls  to  DQXGS 


3 


NSWCDD/TR— 05/90 


IF  (P  <  Pmax)  THEN 
Amax  =  A,  Pmax  =  P 
ENDIF 
GOTO  35 
ELSE 

Amin  =  A,  Pmin  =  P 
GOTO  25 
ENDIF 

25  CONTINUE 
ELSE 

DO  30  J  =  1,  NJ 

A  =  Amax  —  J  *  DaNJ 
CALL  DQ(A,  P,  I) 

IF  (P  <  0)  THEN 

IF  (P  >  Pmin)  THEN 
Amin  =  A,  Pmin  =  P 
ENDIF 
GOTO  35 
ELSE 

Amax  =  A,  Pmax  =  P 
GOTO  30 
ENDIF 

30  CONTINUE 
ENDIF 

35  r  =  DZERO(F2,  Amin,  Amax,  AERR,  RERR)  !DZERO  finds  best  approx  for  rQ. 

The  double-precision  function  DZERO  referred  to  in  the  last  line  is  taken  from  [8,  p.  151]. 
It  finds  a  root  of  the  function  F2  =  P(r)  —  PQ,  where  the  root  rQis  bounded  by  Amin  and 
Amax.  The  absolute  and  relative  errors,  specified  by  the  user,  are  given  by  AERR  and 
RERR.  They  are  presently  set  at  0  and  10”12,  respectively.  Note  that  F2  is  evaluated  by 
a  call  to  DQ,  which  in  turn  calls  DQXGS,  which  calls  the  subprogram  FI  to  evaluate  the 
function  FI  from  (4). 

III.  DETAILS  OF  THE  FORTRAN  77  SUBROUTINE  INVELP 

As  stated  in  the  Introduction,  IELP.FOR  contains  38  subprograms  that  are  used  to  de¬ 
termine  an  acceptable  value  r  for  rQ,  given  P0,h,  k,  u,  v.  All  the  routines  are  designed  as 
double-precision  and  itemized  at  the  end  of  this  section.  The  master  routine,  the  first  listed, 
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has  the  Call:  INVELP(P0,  h,  k,  u,  v,  r,  P,  I).  The  first  5  variables  have  been  defined  in  Section  I 
and  are  input.  The  variable  r  ~  rQ  is  the  desired  output,  with  P  =  P(r)  ~  0.  The  location  I 
contains  the  number  of  calls  to  the  routine  DQ;  it  can  be  as  high  as  55,  provided  (16)  below 
is  satisfied.  The  average  number  of  calls  is  12.  Accuracy  of  the  output  is  discussed  below. 

The  input  constraints  are  : 

10-20  <  PQ  <  1  -  lO-11,  0  <  h  <  108,  0  <  k  <  108,  10“6  <  u,  v  <  106.  (16) 

INVELP  is  designed  to  give  r  to  8  significant  digits  whenever  possible.  However,  because  of 
the  approximately  15-decimal-digit  word  length,  this  will  not  be  achieved  if  H  or  K  is  ex¬ 
tremely  large.4  The  accuracy  in  r  will  generally  be  reduced  for  values  outside  the  inequalities 
in  (16). 

Some  checks  are  carried  out  on  the  input  for  INVELP.  It  is  required  of  the  input  that 
0  <  PQ  <  1.  If  either  of  these  inequalities  is  violated,  r  is  set  to  —2.  Also,  if  either  u  or  v  is 
not  positive,  r  is  set  to  —1.  In  either  case,  with  r  negative,  an  immediate  exit  is  made  from 
INVELP. 

The  third  routine  listed  is  DQ,  which  has  the  Call:  DQ(r,  PX,  II)  where  the  first  variable 
is  input  (with  h,  k,  u,  v  stored  in  common)  and  the  output  is  PX  =  P(r).  DQ  calls  the 
integration  routine  DQXGS  from  [8,  p.  531]  to  evaluate  P(r)  by  (3).  DQXGS  calls  the 
fourth  listed  subprogram,  namely  the  external  function  Fl(s),  which  yields  the  integrand  of 
(3)  for  a  given  argument  s.  In  addition,  the  relative  error  specified  for  DQXGS  is  set  at 
EPSREL  =  5  •  ICC15.  The  possibility  of  reducing  the  integration  interval  in  (3)  is  discussed 
in  Appendix  A  as  carried  out  in  DQ. 

The  fifth  listed  routine  is  GRUB(Rg),  which  gives  an  early  estimate  Rg  ~  rQ.  It  requires 
as  input:  PQ  ,  h,  k,  u,  v,  which  are  accessed  by  a  common  statement.  The  Grubbs’  analysis  for 
this  routine  is  given  in  Appendix  B. 

The  sixth  listed  routine  is  PSQR,  which  was  discussed  in  Section  II.  It  has  the  Call: 
PSQR(Rsq,  Psq,  Rg,  I),  where  Rg,  taken  from  GRUB,  is  an  initial  estimate  for  the  output 
Rsq;  the  other  inputs  are  given  in  common.  The  output  quantity  Rsq  is  the  second  estimate 
for  rQ.  The  other  output  Psq  contains  the  probability  over  P(Rsq)  —  PQ  =  P(Rsq)  (see 
Section  II).  I  contains  the  number  of  (N-R)  iterations  required.  The  analysis  is  discussed  in 
Appendix  C. 

The  remaining  32  routines,  given  in  IELP.FOR,  are  taken  from  [8].  They  are  supporting 
routines  for  the  first  six  routines  discussed  above. 


4 For  example,  if  PD  =  10“10,  h  =  100,  k  =  10s,  u  =  1,  v  =  10-6  =>  rQ  ~  100000000.0000410 
Po  =  .9999999999,  h  =  100,  k  =  10s,  u  =  1,  v  =  10“6  =4>  rD  ~  100000000.0000592 
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1  SUBROUTINE  INVELP(PP,  HI,  Kl,  SI,  S2,  R,  P,  I) 

2  DOUBLE  PRECISION  FUNCTION  F2(X)  !  NEED  ED  FOR  DZERO  OF  INVELP. 

3  SUBROUTINE  DQ(R,PX,II) 

4  DOUBLE  PRECISION  FUNCTION  F1(X)  !FUNCTION  FOR  DQXGS  ROUTINE. 

5  SUBROUTINE  GRUB(R)  ! INITIAL  ESTIMATE  FOR  R  FROM  GRUBBS’  APPROX. 

6  SUBROUTINE  PSQR(R,  P,  Rl,  I)  !1ST  APPROX,  Rl,  USES  GRUBBS’  ESTIMATE  Rg. 

7  INTEGER  FUNCTION  IPMPAR  (I) 

8  DOUBLE  PRECISION  FUNCTION  DPMPAR(I) 

9  DOUBLE  PRECISION  FUNCTION  DEPSLN(L) 

10  DOUBLE  PRECISION  FUNCTION  DXPARG(L) 

11  DOUBLE  PRECISION  FUNCTION  REXP(X) 

12  DOUBLE  PRECISION  FUNCTION  ALNREL(A) 

13  DOUBLE  PRECISION  FUNCTION  RLOG(X) 

14  DOUBLE  PRECISION  FUNCTION  ERF(X) 

15  DOUBLE  PRECISION  FUNCTION  ERFC1(IND,  X) 

16  DOUBLE  PRECISION  FUNCTION  DERF(X) 

17  DOUBLE  PRECISION  FUNCTION  DERFC(X) 

18  DOUBLE  PRECISION  FUNCTION  DERFCO(X) 

19  DOUBLE  PRECISION  FUNCTION  ERFI(P,  Q) 

20  DOUBLE  PRECISION  FUNCTION  DAERF(X,  H) 

21  SUBROUTINE  PNI(P,  Q,  D,  W,  IERR) 

22  DOUBLE  PRECISION  FUNCTION  GAMMA(A) 

23  DOUBLE  PRECISION  FUNCTION  GLOG(X) 

24  DOUBLE  PRECISION  FUNCTION  GAM1(A) 

25  DOUBLE  PRECISION  FUNCTION  GAMLN(A) 

26  DOUBLE  PRECISION  FUNCTION  GAMLNl(A) 

27  SUBROUTINE  GRATIO(A,  X,  ANS,  QANS,  IND) 

28  DOUBLE  PRECISION  FUNCTION  RCOMP(A,  X) 

29  SUBROUTINE  GAMIN V(A,  X,  X0,  P,  Q,  IERR) 

30  DOUBLE  PRECISION  FUNCTION  DZERO(F,  AX,  BX,  AERR,  RERR) 

31  SUBROUTINE  DQPSRT(LIMIT,  LAST,  MAXERR,  ERMAX,  E  ,  IORD,  NRMAX) 

32  SUBROUTINE  DQELG(N,  EPSTAB,  RESULT,  ABSERR,  RES3LA,  NRES, 

*  EPMACH,  OFLOW) 

33  SUBROUTINE  DQXGS  (F,  A,  B,  EPSABS,  EPSREL,  RESULT,  ABSERR,  IER, 

*  LIMIT,  LENIW,  LENW,  LAST,  IWORK,  WORK) 
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34  SUBROUTINE  DQXGSE(F,  A,  B,  EPSABS,  EPSREL,  LIMIT,  RESULT,  ABSERR, 

*  IER,  A,  B,  R,  E,  IORD,  LAST,  VALP,  VALN,  LP,  LN) 

35  SUBROUTINE  DQXCPY(A,  B,  L) 

36  SUBROUTINE  DQXLQM(F,  A,  B,  RESULT,  ABSERR,  RESABS,  RESASC,  VR,  VS, 

*  LR,  LS,  KEY,  EPMACH,  UFLOW,  OFLOW) 

37  SUBROUTINE  DQXRUL(F,  XL,  XU,  Y,  YA,  YM,  KE,  Kl,  FV1,  FV2,  LI,  L2) 

38  SUBROUTINE  DQXRRD(F,  Z,  LZ,  XL,  XU,  R,  S,  LR,  LS) 
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REFINED  LIMITS  OF  INTEGRATION  FOR  P 

The  objective  in  this  appendix  is  to  give  the  analysis  used  in  DQ  to  reduce,  whenever  possible, 
the  interval  of  integration  in  (3).  For  convenience,  (3)  is  reproduced. 

i  rU+K 

P  =  irr  /  F1(=) ds.  (A-1) 

2  J H-R 

where 

Fl(s)  =  exp(— s2)  aerf  [K,  (u/vj^R2  -  (s  -  H)2  ].  (A-2) 

For  P  =  P0,  we  wish  to  find  values  for  a  and  (3,  where 

P  =  P0  =  —^  Fl(s)  ds  +  7r~/=  /  Fl(s)  ds  +  /  Fl(s)ds,  (A-3) 

2  V71  J  H-R  ^  Vn  J  a  2  ^TT  Jp 

and 

H  -  R  <  a  <  (3  <  H  +  R,  (A-4) 

such  that  the  first  and  third  integrals  on  the  right  of  (A-3)  are  negligible  relative  to  the 
second.  Of  course,  in  some  cases  a  =  H  —  R  and  (3  —  H  +  R  so  that  the  integration  limits 
are  left  unchanged. 

Considering  the  first  integral  and  using  the  fact  that  0  <  aerf(a,  b)  <  2,  for  positive  a  and 
b,  one  obtains  for  a  desired  accuracy  specified  by  e 

— -=  j  Fl(s)  ds  <  —=  [  exp(— s2)  ds  <  — ^exp(— a2)  =  PQe.  (A-5) 

2  V  J  H— R  V71  j  H-R  V 7r 

Therefore, 

a  =  a0  =  —  J — log(.5  \/7rP0e/R)  (3  —  (30  —  —  a0.  (A-6) 

Further  efforts  are  made  to  increase  a,  since  F1(H  —  R)  =  F1(H  +  R)  =  0,  by  computing 
a  sequence  of  new  values  cq  =  cq_i  +  jA,  j  =  1,2,  ...,N;  A  =  (/30  —  «o)/N  as  long  as 
Fl(aj_1)  <  PQ  e.  The  procedure  is  continued  until,  for  some  j=J,  the  last  inequality  is 
violated.  Then  the  value  for  a  is  taken  as  aj_i.  Actually  the  algorithm  is  refined  so  that 
two  subintervals  of  A  are  also  used.  A  similar  procedure  attempts  to  improve  the  value  of 
A)  for  /3. 

It  is  worth  noting  that  if  H  —  R  >  —  a0>  then  P  is  set  to  zero.  In  DQ,  at  present,  N  =  50, 
and  e  =  10-11. 
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GRUBBS’  APPROXIMATION  FOR  rQ 

The  approximation  for  rQ  derived  by  Grubbs  [5]  is  often  a  very  good  estimate.  His 
approximation  follows,  with  P  =  PQ  given. 

Let 

a2  =  u2  +  v2  (B-l) 

M  =  1  +  h2/a2  +  k2/cr2  (B-2) 

V  =  2  {  (ir>2)2  +  (v2/ a2)2  +  2  [(u2/ a2)(h2/cr2)  +  (v2/a2)(k2/0-2)  ]  }.  (B-3) 

Grubbs’  analysis  now  calls  on  the  incomplete  gamma  function,5  namely 

P(A,x)  =  fjjT  e_ttA_1dt,  A  >  0,  x  >  0,  (B-4) 

where  V  =  P0,  and 

A  =  M2/V,  (B-5) 

x  =  (M/V)(r2/a2).  (B-6) 

Knowing  A  and  PQ,  the  gamma  inverse  routine  GAMINV  [8,p.  85]  is  used  to  find  x.  Then 
Grubbs’  estimate  for  rQ  is  obtained  by  solving  (B-6)  for  r, 

r  =  \J  x  (V/M)  a2  =  rG.  (B-7) 


5 The  incomplete  gamma  function  is  related  to  the  chi-squared  distribution  [l,p.  262]. 
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APPENDIX  C 

PROBABILITY  OVER  A  SQUARE  CIRCUMSCRIBING  A  CIRCLE 
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PROBABILITY  OVER  A  SQUARE  CIRCUMSCRIBING  A  CIRCLE 

The  objective  in  this  appendix  is  to  first  give  an  expression  for  the  probability,  Psq,  of  a 
shot  falling  under  a  bivariate  uncorrelated  normal  distribution  with  mean  zero  and  standard 
deviations  u  and  v,  in  a  square,  S,  centered  at  (h,  k)  with  sides  parallel  to  the  coordinate 
axes  and  of  length  (2a).  Then  a  brief  description  follows  of  the  algorithm  used  to  find  (a) 
given  Psq,  h,  k,  u,  v. 

The  expression  for  Psq  is  given  by 


Psq  v/2^ 


u 


/»h+a  -i  /»k+a 

/  exp[— (x/u)2/2]  dx  —f== 

lh— a,  V  27 r  V  J k— a 


exp[— (y/v)2/2]  dy.  (C-l) 


Normalizing,  as  was  done  on  page  1  of  the  main  text,  let  x 
becomes 


P  —  _ 

sci  4 


r-H+Au 


r-K+Av 


-4=  /  exp(— s2)  ds  / 

,V7r  Ju-Au  V17  j K-Av 


(\/2u)s,  y  =  (V2)vt,  (C-l) 


exp(— t2)  dt 


(C-2) 


or 

Psq  =  ^  aerf(H,  Au)  aerf(K,  Av),  (C-3) 

where  H  =  h/(y/2  u),  Au  =  a/(y/2  u),  K  =  k/(y/2  v),  Av  =  a/(y/2  v).  The  function  aerf  is 
defined  by  (5);  it  is  evaluated  numerically  by  calling  the  subroutine  DAERF  from  [8,  p.  51] . 

We  start  by  setting  Psq  =  PQ,  and  getting  an  initial  estimate  for  (a)  from  the  Grubbs 
estimate  for  rQ  (see  Appendix  B).  Then  a  halving  procedure  is  invoked,  which  keeps  (a) 
in  the  interval  (amin,  amax),  where  amin  (amax)  yields  a  value  for  Psq  <  (>)  PQ.  When 
amax  —  amin  <  .001  *  amax,  the  Newton-Raphson  iterations  (N-R),  with  an  initial  value  for 
(a)  of  (amax  +  amin)/2,  replaces  the  halving  process. 

The  equation  that  governs  the  (N-R)  is  given  by 


an+i  =  an  -  (Pn  -  P0)/DPn,  n  =  1,  2,  . . . ,  N  ,  (C-4) 

where  ai  =  (amax  +  amin)/2,  Pn  denotes  Psq  from  (C-3)  evaluated  at  a  =  an,  and 
DPn  denotes  the  derivative  of  Psq  ,with  respect  to  (a),  evaluated  at  an,  namely 

DPn  =  —= — exp(H  —  Au)2  [1  +  exp(— 4  H  Au)]  aerf(K,  Av) 

V2  u 

+  —= —  exp(K  —  Av)2  [1  +  exp(— 4K  Av)]  aerf(H,  Au),  a  =  an.  (C-5) 
V2  v 

The  (N-R)  iterations  are  terminated  at  n  =  N  <  60  with  a  =  ax+i,  where 

|ax+i  —  &n|  5;  5 ax+i  10  (C-6) 

It  is  worth  noting,  since  (a)  is  also  the  radius  of  the  circle  circumscribed  by  S,  and  v^a 
is  the  radius  of  the  circle  circumscribing  S,  a  <  rQ  <  \/2a. 
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TABULATION  OF  r  AS  A  FUNCTION  OF  P,  h,  k,  v,  (u=l) 

r  denotes  the  radius  of  a  circle  in  the  xy  plane,  with  center 
at  (h,  k),  that  contains  P%  of  a  cumulative  normal  uncorrelated 
bivariate  distribution  with  mean  zero,  standard  deviations  u  and  v. 


TABLE  D-l.  TABULATION  OF  r  WITH  u  =  1. 

LAYOUT  OF  TABLE 

A  tabulated  value  of  h  is  constant  and  located  at  the  top  of  each  page. 
Tabulated  values  of  k  are  given  along  the  second  row. 
Tabulated  values  of  v  are  specified  down  the  Erst  column. 
Tabulated  values  of  P  are  given  down  the  second  column 

The  ranges  of  the  input  variables  are  (u  =  1): 


h 

=  0.0, 

0.5, 

1, 

2, 

3, 

4, 

5, 

6, 

8, 

10, 

20, 

50, 

120,  500 

V 

=  1, 

2, 

3, 

4, 

5, 

6, 

8, 

10 

k 

=  0.0, 

0.5, 

1, 

2, 

3, 

4, 

5, 

6, 

8, 

10, 

20, 

50, 

120,  500 

p 

=  .01, 

.05, 

.15, 

.30, 

.50, 

.70, 

.90, 

.95, 

.99, 

.999 

Example  1  :  h  =  10.0  v  =  8.0  P  =  0.05  k  =  20.0  r  =  12.1493.  See  page  D-24  for  r. 

Example  2  :  h  =  2.0  v  =  0.5  P  =  0.90  k  =  4.0  r  =? 

The  value  of  v  =  .50,  for  Example  2,  is  not  an  input  value  in  the  table;  nevertheless  r  can 
still  be  found.  By  (3),  P  remains  unchanged  as  long  as  the  ratios  specified  by  R,  H,  K,  u/v 
are  held  constant.  Therefore, 

r  =  ri/ui,  hi  =  4,  kx  =  8,  m  =  2,  v1  =  1. 

Interchanging  hL,  ki  and  Ui ,  y1 : 

r  =  (r2/u2)/2,  h2  =  8,  k2  =  4,  u2  =  1,  v2  =  2. 

The  input  values  for  Example  2,  which  allow  a  table  lookup,  are: 

Example  2  :  h2  =  8.0  v2  =  2  P  =  0.90  k2  =  4.0  r2  =  10.7627, 
with  r  =  r2/2  =  5.3814.  See  page  D-21  for  r2. 
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